
% n = size(mesh.f, 1) * 3;

for i = 1 : 200
%     theta = rand(n, 1);
    pos = rand(3, 1);
    
    diff =  getDeformationGradientFD( pos(1), pos(2), pos(3), mesh, theta ) - evaluateInteriorDeformationGradient( pos(1), pos(2), pos(3), mesh, theta );
    
    check = max(max(diff));
    
    if check > 1e-2
       check 
    end
    
    i
end